

library(tseries)
library(MASS)

# On peut recuperer les indices europeens
data(EuStockMarkets)
cac <- EuStockMarkets[,"CAC"]
dax <- EuStockMarkets[,"DAX"]
smi <- EuStockMarkets[,"SMI"]
ftse <- EuStockMarkets[,"FTSE"]

# SP500 contient les log returns

# une fonction idiote pour calculer les log-returns 
returns = function(X) {
  return (log(X / lag(X)))
}

X <- returns(cac)
layout(matrix(1:2, 2, 1))
plot.ts(cac, main="index")
plot.ts(X, main="log-returns")

X <- SP500
layout(matrix(1:2, 2, 1))
plot.ts(exp(cumsum(X / 250)), main="index")
plot.ts(X, main="log-returns")

# ACF des log-returns et du carre des log-returns
layout(matrix(1:2, 2, 1))
acf(X, ci=0, main="ACF de X")
acf(X^2, ci=0, main="ACF de X^2")

# analyse univariee
layout(matrix(1:2, 2, 1))
plot(density(X), main="Estimateur de la densite des log-returns")
qqnorm(X, main="QQ Plot")
qqline(X)

# on fitte un GARCH(1, 1)
fit <- garch(X - median(X), order=c(3, 2))
summary(fit)
plot(fit)

layout(matrix(1:2, 2, 1))
ts.plot(X, main="returns")
ts.plot(fit$fitted.values[,1], main="Estimates of the conditional variance")

## etude des residus
res = fit$residuals[!is.na(fit$residuals)]

jarque.bera.test(res) # gausiannite regetee fortement...
Box.test(res, type="Box-Pierce")  # test d'independence
Box.test(res, type="Ljung-Box")  # test d'independe

layout(matrix(1:2, 2, 1))
acf(X^2, main="returns^2")
acf(res^2, main="residuals^2")

Box.test(rnorm(100), type="Ljung-Box")

denres = density(res)
plot(denres)
lines( denres$x, dnorm(denres$x, mean(res), sd(res)), col="blue")

qqnorm(res)
qqline(res)

## Illustration de l'effet de levier
Xplus <- na.remove( (X[2:length(X)])[ sign(X)==1 ] )
Xminus <- na.remove( (X[2:length(X)])[ sign(X)==-1 ] )

layout(matrix(1:2, 2, 1))
plot.density(density(Xplus))
plot.density(density(Xminus))

qqplot(Xplus, Xminus)

# on dessine la qqline a la main
q1 = as.numeric(quantile(Xplus, probs=c(0.25, 0.75)))
q2 = as.numeric(quantile(Xminus, probs=c(0.25, 0.75)))
a <- diff(q2) / diff(q1)
b <- q2[1] - a * q1[1]
x = seq(min(sort(Xplus)), max(sort(Xplus)), 0.1)
lines(x, a*x + b, col="red")

qqnorm(Xplus)
qqline(Xplus)



# prediction

condvar <- fit$fitted.values[, 1]

prediction <- predict(fit)[,1]

layout(matrix(1:2, 2, 1))
plot.ts(prediction)
plot.ts(condvar)

## predict calcule juste la var conditionnelle
## si on veut vraiment predire, il faut utiliser
## l'option genuine=TRUE

# on enleve les 100 dernieres observations
b = 100

training = X[1:(length(X)-b)]
test = X[(length(X)-b+1):length(X)]

fit <- garch(training, order=c(1, 1))
condvar1 <- predict(fit, test)

fit2 <- garch(X, order=c(1, 1))
condvar2 <- predict(fit2)


# layout(matrix(1:2, 2, 1))
# pour illuster la difference entre prediction et filtrage
plot(condvar1[,1], type="l")
lines(condvar2[(length(X)-b):length(X),1], col="blue")

# prediction de la vol du SP500
fit <- garch(X - median(X), order=c(1, 1))


condvar2 <- predict(fit, test, genuine=TRUE)



pred <- predict(fit, genuine=FALSE)[,1]
predg <- predict(fit, genuine=TRUE)[,1]


b = 100
last <- pred[ (length(pred)-b):length(pred)]
plot.ts(last)
lastpred <- (length(last)-1):length(last)
lines( lastpred, last[lastpred], col="red")




